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SOLVING UPWIND-BIASED DISCRETIZATIONS: DEFECT-CORRECTION 

ITERATIONS 


BORIS DISKIN* AND JAMES L. THOMAS* 

Abstract. This paper considers defect-correction solvers for a second order upwind-biased discretization 
of the 2D convection equation. The following important features are reported 

1. The asymptotic convergence rate is about 0.5 per defect-correction iteration. 

2. If the operators involved in defect-correction iterations have different approximation order, then 
the initial convergence rates may be very slow. The number of iterations required to get into the 
asymptotic convergence regime might grow on fine grids as a negative power of h. In the case of a 
second order target operator and a first order driver operator, this number of iterations is roughly 
proportional to A -1 / 3 . 

3. If both the operators have the second approximation order, the defect-correction solver demonstrates 
the asymptotic convergence rate after three iterations at most. The same three iterations arc required 
to converge algebraic error below the truncation error level. 

A novel comprehensive half-space Fourier mode analysis (which, by the way, can take into account the 
influence of discretized outflow boundary conditions as well) for the defect-correction method is developed. 
This analysis explains many phenomena observed in solving non-elliptic equations and provides a close 
prediction of the actual solution behavior. It predicts the convergence rate for each iteration and the 
asymptotic convergence rate. As a result of this analysis, a new very efficient adaptive multigrid algorithm 
solving the discrete problem to within a given accuracy is proposed. Numerical simulations confirm the 
accuracy of the analysis and the efficiency of the proposed algorithm. The results of the numerical tests are 
reported. 

Key words, convection, upwind-biased discretization, defect- correction iteration 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. This is the first in a series of papers analyzing the efficiency of different iterative 
algorithms solving upwind- biased discretizations of the convection operator. The model problem we study 
in this paper is the 2D constant coefficient convection equation 


( 1 - 1 ) LU=(a-v)u = F(x,y), 

where a — ( 01 , 0 , 2 ) is a given vector. 

The solution U(x,y) is a differentiable function defined on the unit square (x, y) £ [0,1] x [0,1]. In 
this paper, we deal mostly with the homogeneous equation F(x,y) = 0. Exceptions when non-homogeneous 
right-hand side functions F(x i y) are considered will be emphasized specifically. 
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Let (f> be the non-alignment angle (another name which is common in CFD is the angle of attack), i.e., the 
angle between the vector a and the positive direction of the x axis; t = tan <f> = a 2 /ai is the non-alignment 
parameter. For simplicity we assume dj > a 2 > 0 and, therefore, 1 > t > 0. Then, Eq. (1.1) can be rewritten 
as 

( 1 - 2 ) d^U = F(x,y)/y/al + aZ, 

where f = is a variable along the characteristic of Eq. (1.1). 

Eq. (1.1) is supplied with Dirichlet boundary conditions at the inflow boundary x = 0 and periodic 
conditions in the y direction. 

( L3 ) U{®,y) = 9{y), U(x,y) = U(x,y + 1), 

where g(y) is a given function. 

This problem is discretized on the 2D Cartesian uniform grid with meshsize h in both the x and y 
directions. Let u lui2 be a discrete approximation to the solution U(x,y) at the point (x, y) = (iih,i 2 h). 
Then, the second order accurate discretization corresponding to the Van Leer’s scheme with k = 0 is defined 
as 


(1.4) 


L h u x 


\+t 2 h ^ ^ w *i + l,t2 “b 3w ii,t2 3u ii — l,t 2 d" u ij -2,i 2 ^ 

2v /l ( ( 3u N,i 2 “ 4 u W 1,*2 + u N-2 ,i 2 ) 

+t(3u Nt i 2 - 4ujv,i a _i + Ujv,i 2 -2)^ = /Af,i 2 ; 
h = 1,2,...AT-1, i 2 = 1,2,. ..AT, N=l/h- 
u o,i 2 = 9 {hh), u-i,i 2 = g'ihh). 


L h u N , i2 = 


In computing L h at nodes with z 2 = 1, z 2 = 2, and 12 = N the vertical periodicity is employed. 
The outflow boundary conditions at i\ = TV are discretized by the second order upwind scheme. The 
discretization of the right-hand side function is f iui2 = F(Uh, i 2 h)/y/af+ a\. Function g'(y) is an additional 
numerical boundary condition. In model problems where the exact solution U(x, y) is known, one can define 
9*{y) “ U(—h, y). 

The discrete scheme (1.4) is upwind biased but not a pure upstream scheme since for defining the 
operator value at the point ( 11 , 12 ), the solution values at the downstream points (i x + 1 , 12 ) and + 1) 

are required. 


1.1. Defect- Correction Schemes. The subject of this paper is the defect- correction method (see [6]) 
which is widely used to iterate non-upwind discretizations. Let our target discrete problem be 


(1-5) L h u llM =f iui2 , 

where L h is a upwind- biased discretization of the convection operator (e.g., (1.4)). 

Let Lg t be a stable (say, the first or the second order pure upwind) discretization of the same convection 
operator and Ui u i 2 be the current solution approximation. Then the improved approximation u iui2 is 
calculated in the following two steps. 
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1. The correction Vi lt i 2 is calculated by marching operator L^ t with a right hand side being represented 
by the residual computed for the target operator L h with the current approximation. The inflow 
boundary conditions for v are initialized with the zero values. 

( i * 6 ) ^ ^ 0,»2 

2. The current approximation is corrected 

(T7) = ^u>i2 

The operator L^ t is called the driver operator. The steps (1.6)- (1.7) can be repeated until the desired 
accuracy is reached. 

In many papers (e.g., [5]), authors studying the defect-correction iteration for non-clliptic problems ob- 
served a slow convergence or even divergence in some common error norms for the initial iterations and 
good asymptotic convergence rates afterwards. This behavior is quite different from that demonstrated 
by the defect-correction method in solving elliptic problems where the asymptotic convergence rate is the 
slowest one. We realized that this “non-elliptic” feature is explained by some properties associated with the 
cross- characteristic interaction (e.g, dissipation) in the operators involved in the defect- correction iterations. 
Namely, this cross-characteristic interaction and the frequency of an incoming component define the pene- 
tration distance (also termed “survival distance” in [4]) of this component. The penetration distance is the 
distance from the inflow boundary on which the discrete solution of the homogeneous problem reasonably 
approximates the continuous one (i.e., the norm of the relative discretization error is essentially less than 
1). The ratio of penetration distances of the operators L h and L^ t determines the number of defect- correction 
sweeps required to get into the asymptotic convergence regime. Moreover, comparison of the penetration 
distances is a constructive and convenient tool allowing one to decide which grid is appropriate for a given 
problem, provided information is known regarding flow direction variations, the frequency content of inflow 
boundary conditions, the characteristic size of the domain, and the accuracy desired. 

The analysis of the first differential approximations (FDA) (see [7]) to the operators involved in the 
defect- correction iterations is presented in Sec. 1.2. It provides us with some qualitative description of 
corresponding penetration distances. This analysis concludes that if the operators L h and L^ t have different 
approximation orders (the case described in Sec. 2) then the efficiency of defect-correction iterations is, 
actually, grid dependent. It means that the maximal number of sweeps which might be required to reach 
the asymptotic convergence rate (or to reduce the algebraic error to the truncation error level) on fine grids 
is larger than on coarse grids. This phenomenon relates to the fact that some smooth error components are 
poorly approximated by the driver operator. In other cases, when the operators L h and L% t have the same 
approximation order (see Sec. 3), the defect- correction method can serve as a very efficient solver. 

1.2. Poor Characteristic Component Approximation. The non-ellipticity of the operators to be 
analyzed introduces a new issue in the standard discretization analysis common for elliptic operators. The 
discretization error of a discrete elliptic operator regarding a given component is defined by (1) the operator’s 
approximation order and (2) frequencies of the component. For non-elliptic operators, the main factor is 
(3) the penetration distance d of this component. The penetration distance of an incoming component 
depends on its frequency cj, i.e., d — d(u>). Note, that in the homogeneous differential problem (1.1)- 
(1.3) all the incoming oscillations are translated along the characteristics without changing their phases 
and amplitudes. However, on a grid non-aligned with the characteristic, any discretization unavoidably 
introduces some numerical cross-characteristic interaction which damps amplitudes and/or shifts phases of 


3 



the incoming components. A quantitative measure of this numerical cross-characteristic interaction is the 
coefficient of the largest cross-characteristic derivative appearing in the first differential approximation to the 
operator under consideration. When the cross-characteristic interaction is strong, there are no approximation 
properties provided by the discrete operator beyond some o(d) neighborhood of the inflow boundary. A 
discretization is considered to be accurate with respect to a given incoming component if and only if the 
penetration distance of this component is comparable with (or exceeds) the characteristic size of the domain. 
Of course, the penetration distance of a given component increases on grids with smaller mesh sizes. The 
coarsest grid on which the penetration distance approaches the characteristic size of the domain can be 
considered as the optimal grid for resolving this component. The expected discretization behavior (i.e, the 
discretization error of a p-th order accurate discrete operator is reduced by factor 2 P when the grid mesh size 
is refined to h/2) is observed only for grids which provide an accurate resolution for all the essential inflow 
components. 

Let L h be an accurate discretization with respect to some particular inco min g component. To approxi- 
mate the solution of L h operator by solving some less accurate operator L Bt (with a correspondingly shorter 
penetration distance), one has to iterate L s t as many times as needed to attain accuracy up to the L h 
penetration distance. The following analysis shows that the required number of iterations depends on some 
power of N , where N is the number of grid points in the characteristic direction. To be precise, the number 
of iterations is proportional to N where p and r are the approximation orders of operators L h and L Bt 
respectively. Below, we derive the predicted dependence for a particular case where p — 2 and r = 1 . 

Before proceeding with the analysis, let us introduce a definition. The components which are much 
more smooth in the characteristic direction than in other directions are referred to as the characteristic 
components. 

The target operator L h approximates the differential operator L from ( 1 . 1 ) with second order accuracy. 
It means that its first differential approximation (FDA) is 

FDA(L h ) = - C 2 h> (d vm + d c B 2 (dt,d v )), 

where C 2 is a constant, 77 = is the cross-characteristic coordinate and is a linear combina- 

tion of the second-order derivatives in £ and r?. For characteristic components (in terms of which » <%), 
this approximation is simplified to 

(1-8) FDA(L h ) ndt-C 2 h 2 d VT)v . 

Let the driver operator L^ t have the first order accuracy. Then, its FDA taken for the characteristic 
components is 

( L9 ) F£>a(l*) W0 c -C,MU, 

where C\ is a positive constant. 

We are going to perform a half-space mode analysis for the first differential approximations to the 
operators L h and L^ t and demonstrate in this way that the driver operator (and the whole defect- 
correction iterative process) is likely to approximate poorly some smooth characteristic components of the 
second order solution. 

Following [ 2 ] and [4], the half-space mode analysis presented in this section is focused on approximating 
the characteristic components. It considers the discretizations of the homogeneous equation ( 1 . 1 ) on the half 
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space |(x, y) : x > Oj with boundary conditions (on x = 0) being represented by one Fourier mode e tujy at 
a time. The purpose of this analysis is to estimate the penetration distance as a function of the incoming 
oscillations. 

For the driver operator, we seek a bounded differentiable function 0(x, y) satisfying the following equation 
and boundary conditions 


(1.10) drf - Cxhdmt = 0; <t>(0 ,y) = e iu}y . 

The exact solution can be written out as 

<j) — e Cih0 2 Z+Pri, 


where (3 = a + ib is a complex number with a and b satisfying the system of the algebraic equations 

I C\h(a 2 — b 2 )t -f a — 0; 

( ‘2C\h(ibt -f- b — \/TT t^uj. 

From the system, a = O(h) and b = yj 1 + t 2 uj + 0(h ). Therefore, the leading term of the bounded solution 
to (1.10) is 

(1,11) (f) ^ e ~C 1 h(l+f 2 )u; 2 £+tv / l+« 2 u>T^ 


The factor e Vi+t 2 {iu>t)) d escr ibes the exact solution of the continuous problem while the factor e -CiMi+* 2 )^ 2 £ 
is the influence of the numerical cross-characteristic interaction. In the case of the first order driver, this 
is a dissipation which damps the amplitude. From (1.11), one can observe that the penetration distance 
on which this damping becomes 0(1) is proportional to . In a similar way one can derive 

the penetration distance of a second order scheme, which is proportional to . For even order 

schemes, the only difference is that the numerical cross-characteristic interaction usually affects the phase 
of the incoming component rather than the amplitude. The comparison of the penetration distances d ^ 
and d ^ implies that the number N sweeps of defect-correction sweeps (1.6)- (1.7) needed to approximate the 
incoming component e lu>y to the second order accuracy is 


( 1 . 12 ) 


AT., 


dW 


1 

u)h 


It is obvious that this number increases when uj tends to zero. On the other hand, if a ; is sufficiently small 
(e.g., the component is nearly constant) then the penetration distance for any scheme covers the entire 
domain and the desirable accuracy is achieved. This consideration implies that the worst case is the case 
d W « R and d ^ R , where R is a characteristic size of the domain ( R = \/l + t 2 in our problem). Thus, 
on the given grid, the worst component is uj ~ R~*h~i , and for this component 


N, 


sweeps 



This meshsize dependence was first mentioned in [3]. It is really not very harmful and in many practical 
calculations it can hardly be noticed. However, an accurate choice of data in the numerical experiments 
allows us to observe this behavior. 

Further, in Sec. 2, a defect-correction method with the first order driver is considered. We introduce 
another (nearly) precise discrete half-space mode analysis in Sec. 2.1 which then will be used to predict 
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the accuracy and the efficiency of the defect-correction iterations (tested in Sec. 2.3). Some analytical 
predictions about the asymptotic convergence rate of residual norms are made in Sec. 2.2 and are validated 
in Sec. 2.3. The comparison of the analytical and the numerical results for the defect-correction method 
with a second order driver is presented in Sec. 3. A very efficient adaptive multigrid algorithm yielding the 
approximate solution to a desired accuracy is proposed and tested in Sec. 4. 


2. Defect-Correction Method with First Order Driver (DC1). The driver used in the DC1 

algorithm is the first order upwind discretization of the convection operator (1.1). 


( 2 . 1 ) 




*i = 1,2,. ..JV, 
u o,i 2 = g(hh). 


*2 = 1,2,. ..AT; 


This scheme is stable for downstream marching. In our case of y-directional periodicity, the marching of this 
scheme requires an implicit line-by-line rather than simple pointwise passage. The entire defect-correction 
iterative process has already been defined above in Sec. 1.1. 

2.1. Half-Space Mode Analysis. In this section we exhibit a discrete half-space mode analysis of 
the defect-correction method which is distinct from the FDA analysis presented in Sec. 1.2: it considers 
the discretizations themselves rather than their differential approximation. This tool is much more accurate 
(and cumbersome at the same time). It can be used to explain in detail many phenomena observed in solving 
non-elliptic equations and provides a close prediction of the actual solution behavior. 

This analysis considers each discretization on the half space as it is, while the boundary condition is 
represented by a Fourier component. In this way the original multidimensional problem is translated into a 
ID discrete problem, where the frequencies of the boundary Fourier component are considered as parameters. 
To regularize the half-space problem, the solution is not allowed to grow faster than a polynomial function. 
The goal of this analysis is the comparison with each other of (1) the exact solution of the differential 
problem, (2) the exact solution of the discrete problem and (3) the approximate solutions at different stages 
of the solver. 

2.1.1. Exact Solutions and Discretization Error. Choosing the domain to be {( x,y ) : x > 0}, for 
each Fourier frequency u> 2 , the differential problem (1.2), (1.3) can be reformulated in the following way: 
find function U(x,y) such that 


d^U = it, U{0,y) = e iui2y , 

where = (u q + tu 2 )/ Vl + 1 2 is the characteristic frequency (/% « 0 for characteristic components). The 
exact solution of this problem is U{x,y) = 

The discrete counterpart is 


(2.2) L h u ilth = i/J {e *(«m+n 2 i 2 ), ^ = n _ lA = e i(-n 1+ n 2 i 2 ) , 

where L h is the target discrete operator, Sli = uq/i and D 2 = w 2 /i are normalized frequencies. 
We seek a solution of the discrete problem in the form 
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(2.3) u ilM = 4> tl e in ^. 

Then, the problem (2.2) can be reformulated for <f) tl as 

(2.4) a - 2 (^ 2 ) 0m - 2 + + ao (^2 + ai ^l 2 )0u+i = V^l -f 

(2.5) 0o = l, 0-i=e~ in \ 

where 

a- 2 (« 2 ) = 
a-i(n 2 ) = -j; 

(2.6) a 0 (fi 2 )= | + |^e- i2n2 -5e- in2 +3 + e <n ^; 

ai (f) 2 ) = i. 

The solution to (2.4), (2.5) is given by 

(2.7) 4> h = W 0 e* niil + (l - Wo) (Corj* + Cir* 1 ), 
where ro and r\ are the roots of the cubic equation 


a-2 


“f~ 0—1 (^2^) r ~b do ~b a l 


]r 3 =0, 


satisfying to |r| < 1, 


W 0 = 


Co 


a- 2 (^ 2 ^ 
ro (n - e 


0 — i2n 2 | 
ifli 


y/lT t 2 hi(3 $ 

a^i(ft2^e~ ini -ha 0 ^ 2 ) -bai^ 2 je‘ ni 


0 


Cl = 


gifi 1 

T\ 


{r\ - r 0 ) 
(fo - e ' Ql ) 
2l (r 0 -r 1 ) 


We avoid here considering in details the exceptional cases where either the denominator in the expression 
for Wo turns out to be zero or ro =r\. In these cases, the form of the solution (2.7) remains the same while 
W 0 and/or Cj (j = 0, 1) might turn to some linear functions of i\. 

Thus, the discretization error is calculated as 


^ii/i,i 2 /i) - 

u iiyi 2 ~ 

e ini<1 - <Ai 1 

(l - W 0 ) 

e iniil - Cor* 1 - c ir \ l 


gtf22*2 

e tH 2 i2 . 


7 



2.1.2. DC1 Iteration. Let the boundary conditions be represented by a discrete Fourier mode e in 2 ' 2 . 
Then, the first order driver operator can be rewritten as 




hy/l + t 2 




0 in 2 i 2 


where Ui u i 2 = c^e^ 2 * 2 is an approximate solution and 

d - 1 ^ 2 ) = —1; = 1 + . 

The general form of components appearing at any stage of the defect- correction iteration is 

(2.8) P(n)g a V n2 S 


where P(z*i) is a complex- coefficient polynomial of i\ and q (|(f| < 1) is the base of the given component. This 
form is invariant under all the transformations occurring in the DC1 iteration. In fact, for any component, the 
only part to be changed is the polynomial P(i\ ) which will be referred to as the amplitude of the corresponding 
component q il e 1 ^ 2 12 . This allows us to analyze separately any building block of this algorithm, such as the 
calculation of residual or solving the driver equation. The underlying idea is to use computer capabilities 
already at the step of deriving an analytic representation for the current solution approximation. We actually 
analyze the response of each building block to an input component in the form (2.8). The output of the 
block is formulated in the same form (2.8), except that the block may produce several output components, 
differing in their bases q and/or frequencies ^ 2 * 

Usually, the initial solution approximation satisfying the boundary conditions cannot be represented as 
a sum of a finite number of components in the form (2.8), e.g., the zero approximation inside with given 
(non-zero) boundary conditions at i\ = 0 and i\ = — 1 . In this case, the collection of analytical components 
(2.8) well describes the distant behavior of the approximation while an adjustment is still needed in the 
neighborhood of the boundary. The approximation in the neighborhood is given by an additional pointwise 
component 


f B h e iQ 0 < h < N 0 
[ 0 otherwise, 

where B h is a complex- valued vector of the length N 0 . iV 0 = 0 at the beginning for many reasonable 
initial approximations including (1) zero approximation, (2) solution of the driver equation, and (3) solution 
interpolated from the coarse grid in the framework of a 2-level multigrid solver. Then, Nq is increased by 
1 in each DC1 iteration. The segment 0 < i\ < No we will call the pointwise representation region and the 
vector B will be called the pointwise amplitude while the domain No < i\ will be referred to as the analytic 
representation region. 

Thus, the first necessary step in analyzing how the DC1 iteration affects the given approximation is to 
separate all the components (including the pointwise component) in the approximation. 

Let = P(ii)q l1 e 1 ^ 2 * 2 be a particular component of the current solution approximation. We are 

going to trace all the changes happening with the amplitude P(i x ) of this component in a DC 1 iteration. The 
changes in the pointwise component will be emphasized as well. An example of the DC1 iteration analysis 
will be presented below in this section. 
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Taking Residual. The residual amplitude /?(ii) of the component Vi u i 2 is calculated as 
R(ii) = A - — L=j a_ 2 (i7 2 ) 9 - 2 P(i 1 - 2) + a_i - 1) 

+ao(fi 2 ).P(ji) + ai + 1) , 

where A = i/3 . f if qr ll e tr * 2 * 2 is the right-hand side component (q = e lQl ), otherwise A = 0. 

Let JVJ = Ah + 2 and B lx = 0,ii > No- Then, the pointwise residual function is computed in the 
following way 

a-2 i ^ r — i d" l d~ a o ^^ 2 ^ d - a i (^ 2 ) u+i > 1 5 

a_2 f S 0 + a~i d- ao d- o>i H — 2; 

a-2(n 2 )flii-2 + a ~i + a\ 3 < i\ < N q; 

0, otherwise. 

where 5_i and So are the boundary conditions (at i\ = — 1 and i\ — 0 respectively) in the problem associated 
with the pointwise representation region. 

Correction Calculation.. The amplitude C(i\) of the correction to the component V{ x ^ 2 is calculated 
from the equation 

d_ 1 (n 2 ) 9 “ 1 C(i 1 - 1) + do (n 2 )c(ii) = hy/l+t?R{h). 

If Vi 1 .i 2 is not an eigencomponent for the driver operator (g / di = — do(fi 2 )/d_i(Sl 2 )), then the power of 
the polynomial C(i i) is the same as of the polynomial R(i i); otherwise the C(i\) power is higher. 

To satisfy the zero conditions at the incoming boundary i\ = 0 which accompany the correction equation, 
one has to complete the correction with the eigencomponent Dod^ 1 e lfl2t2 with the amplitude 

D 0 = ~C( 0). 

The pointwise correction Cf* is calculated from the following system of linear equations. 

c% = 0 
iV 0 

d— i (n 2 )cf i t +d 0 (n 2 )c > [; t +1 = hVT+t*R^ +1 ,o < *, < jv'. 

The amplitude D\ of the accompanying eigencomponent is computed by 



New Amplitude.. The new amplitude P{i\) of the component is calculated as 

P(»i)=P(ii) + C , (< 1 ). 

The corrected pointwise amplitude Bi 1 and the new boundary Nq of the pointwise representation region are 

B n = Bi, + CP‘; N 0 = N 0 + l. 

The amplitude D{i\) of the eigencomponent is also changed to D{i\) 

D{i\) — D(ii) d - D 0 d- D\. 
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2.1*3. Discretization of Outflow Boundary Conditions. The discretization of the outflow bound- 
ary conditions can be taken into account as well. Discretized outflow boundary conditions usually imply 
some special discretization stencil different from that in the interior. In order to incorporate this feature 
into the analysis we have to introduce another pointwise representation zone near the outflow boundary In 
other words, an additional pointwise component is required to make the half-space analysis be precise for all 
the y-pcriodic problems. In the analysis tests below (Sec. 2.3.1) we did not implement this adjustments. It 
was realized that the analysis is actually precise even without this special care about the outflow boundary 
conditions. 

2.1.4. Example of Analysis. Let us consider the homogeneous problem (/% = 0) with boundary 
conditions given by (2.5) and the zero initial approximation inside of the half-space 0 < i x . At the beginning 
the only component involved in calculation is the zero-length (N 0 = 0) pointwise component. The problem 
associated with this component is 


0 - 2 ^ 2 ^ tl _ 2 + a_i j — i. + + a\ (n a )B il + i = 0; 

S 0 = 1, S_! =e~ ini 


1. The extended boundary of the pointwise approximation region is Nq = 2. 

2. The residual function is 

R 1 1 ~ i+t 5 (<*-2(^2) 5 - 1 +a_i(ft 2 )s 0 ^; 

^2 = -^ 7 TTF a - 2 ( ih ) So ’ 

„ Rtf — 0, 2 < 

3. A non-zero correction proves to be only at point i\ = 1 

cpt = hV 1 + 

d - 1 (fl 2 ) 

B\ = Cf; 

A^o = i; 

The eigencomponent is introduced with the amplitude 

hy/l + t 2 R^ - do ( fia'jcf 

D x - — ^ . 

d - 1 (^ 2 j 

4. The approximation obtained at the end of the first DC1 iteration is 


D\d\ + B\, 

Did!? , 


Q(i 1) = 


ii = 1; 

»i > 1; 


5. On the next iteration the initial approximation contains the eigencomponent with the amplitude 
D{ii) = D 1 . Therefore, the boundary conditions in the problem associated with the pointwise 
component is changed to 


S 0 = 1 — £>(0), S_! = - D(—l)/di. 
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The calculations can be continued. Approximations obtained in these iterations can always be repre- 
sented as 


i ,i-2 


Q(ix)e in ^\ 


Q(i i) = 


D{ii)d l i -f Bi 1 , 1 < u < N 0 ; 
, A^o < i\\ 


2.1.5* Algebraic and Total Errors. Let additional component Wi u i 2 = W (ii)e tQl11 e lU2t2 be involved 
into the calculation due to the right-hand side function; then, after any iteration the current approximation 
representation is a sum of this component, the cigencomponent d ily { 2 = D^i)^ 1 e tQ2t2 , and the pointwise 
component. Thus, the current solution approximation is 


(2.9) 


^*1 ,12 


= Q(i 


Q(h) 


w(*i)c <niil + Dih)^ +B tl , 
W(h)e m ^ + D{h)^ , 


1 < *i < N 0 - 

N 0 < ii- 


The algebraic error function which is the difference between the exact and approximate solutions of the 
discrete problem is given by 


( 2 . 10 ) 


i in 


Q(ii) ~ (W 0 e in ^ +Cor£ +C 1 rj 1 ) 


1^212 


The total error function defined as the difference between an approximate solution of the discrete problem 
and the exact solution of the differential problem is calculated as 


( 2 . 11 ) 


- U(iih,i 2 h) 


Q(h)-e in 


gifl2*2 


Involving other components, either due to a non- zero initial approximation or since the algorithm itself 
produces additional components, extends the number of items in (2.9) with straightforward changes in the 
algebraic and total error expressions (2.10) and (2.11). 

2.2. Convergence in Residual Norms. The matrix analysis described in this section can be con- 
sidered as the asymptotic case of the half-space analysis (with the discretized outflow boundary conditions) 
exhibited in Sec. 2.1. In this asymptotic regime, the pointwise representation region covers all the domain. 
Let the problem (1.1)-(1.3) be defined on a layer {x,y) € [0,1] x (— oc, +oo) with the input data (functions 
F(x,y) and g(x,y)) such that the function U(x,y) = ${x)e %uiy is the exact solution of the problem. The 
approximate solution of the corresponding discrete problem (1.4) is sought in the form (2.3). The problem 
for the discrete function 0^ is derived similar to the problem (2.4)- (2.5) 


0o = 4>(0), 0- 1 = $(— h), 

(2.12) a- 2 ^2^0ii-2 T G _ 1 ^2^011-1 + Go ^2^0ii + Gi ^2^0»i + l 
= hyj 1 + t 2 ^4> x (iih) + j = — 1; 

6_2^2^0;v-2 4- 6-i ^O2^0 at-i + bo (n 2 )0^ = hy/l + t 2 ^3> x (l) + ftu;<I>(l)^. 
Functions aj ( j = —2, —1,0, 1) are defined in (2.6), 

f>- 2 (n 2 )=4; 
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( 2 . 13 ) 


6-i (^ 2 ) = -2; 

bo(n 2 )= ^ + ^(e-* 2 ^-4e-^+ 3 | 


Let the TV-dimensional vector 4> = ), (zi = 1, . . . , TV) be an approximate solution to (2.12) with the 

algebraic error e = ^e tl j , i.e., 

4 > = 4 > ~ e.! , 

where </> is the exact solution to (2.12). 

The correction J is calculated from the linear system of N equations 

D5 = — Te, 

where Af-by-N banded matrices T and D correspond to the target operator and to the driver operator 
respectively. 


T = 


ao 

a x 

0 

0 

0 


0 

0 

0 

0 

a_i 

ao 

ai 

0 

0 


0 

0 

0 

0 

0-2 

a- 1 

C lQ 

ai 

0 


0 

0 

0 

0 

0 

a-2 

a~i 

ao 

ai 


0 

0 

0 

0 

0 

0 

0 

0 

0 

• • a- 2 

0-1 

ao 

Ol 

0 

0 

0 

0 

0 

... 

0 

6—2 

6-i 

6o 

do 

0 

0 


0 

° ) 





d- 1 

do 

0 


0 

0 





0 

d-i 

do 


0 

0 





0 

0 

0 


d- 1 

do / 






D 


Then, the amplification matrix G e of the defect correction iteration 

-new = 

becomes 

G e = I - D~ l T . 

Since the residual is usually used in practice to monitor the error, we can modify G e to measure the residual 
reduction, as 


G r = TG e T 


-l 


so that 


f new = q f 
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where 


r — —Te 

represents the discrete residual. If we wish to bound the amplification of the residual, we could use the 
spectral radius of G r ( p{G r )) or the Z/2-norm of the matrix G r (||G r || 2 = yj p(G*G r )). The spectral radius 
p{G r ) is usually referred to as the asymptotic convergence rate, i.e., the rate corresponding to a large number 
of iterations. The L2-norm (||G>||2) indicates the “worst” possible convergence rate. 

2.3. Numerical Tests. 

2.3.1. Verification of Analytical Predictions: Discretization Accuracy Test. The discretiza- 
tion accuracy of an operator on the grid of meshsize h for a characteristic component ( w herc 

u)\ + tuj2 ~ 0) is defined by the penetration distance of this component. (Sec Sec. 1.2.) The first differential 
approximation to the target second order accurate discrete operator ( 1 . 4 ) is defined by 

FDA^L ^ = d x 4- td y — -^-j==== 4* tdyyy'j ■ 

For characteristic components, it turns to 

i 2 

( 2 . 14 ) FDA(L h ) =%- (_t 3 + ^ w; , 

12(l + t 2 J 

where the non-alignment parameter t and the characteristic variables f and 77 are defined in Sec. 1 . Notice, 
when t « 1, the coefficient of the third derivative with respect to 77 vanishes in the FDA and the next term 
(the fourth derivative) becomes important. 

In the general case 0 < t < 1, the discretization error of L h ( DE(L h )) can be approximately calculated 
from the asymptotic solution of the half-space problem associated with ( 2 . 14 ). The derived estimate is 

DE(L h ) = e <(u,1I+W23/ ) (l - e~£ 


where 

( 2 . 15 ) 



Definition 1 . We say that a discretization has the accuracy e for the given characteristic component 
urith the inflow y -directional frequency u )2 on the distance S (measured along the characteristic) from the 
boundary if the following inequality holds 


( 2 . 16 ) 


1 — e ^ < e, for £ < 


where d is the normalized penetration distance of the given characteristic component. Note, that e defines 
the relative accuracy, hence 0 < e < 1 and e % 1 indicates very poor accuracy. For operator ( 1 . 4 ), d — ^2 
and the penetration distance of the e-accuracy is estimated from ( 2 . 16 ) as 


( 2 . 17 ) 


82 = |d 2 | arccos(l — e 2 /2). 
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Table 2.1 

Penetration distances (m meshsizes) of the 1% accuracy. 


t 

U)2 

Penetration distances 

target operator 

driver operator 

<si 2) 


s {2] 

°N 


W 

IF* 

°N 

0.2 

2n 

256 

256 

256 

139.034 

139 

139 

4n 

256 

256 

256 

34.758 

34 

34 

8tt 

256 

256 

256 

8.69 

8 

8 

107T 

256 

256 

256 

5.561 

5 

5 

16 tt 

82.564 

81 

81 

2.172 

2 

2 

0.6 

27 r 

256 

256 

256 

34.758 

34 

34 

47T 

256 

256 

256 

8.69 

8 

8 

87T 

256 

256 

256 

2.172 

2 

2 

107T 

169.092 

162 

162 

1.39 

1 

1 

I67 r 

41.282 

37 

37 

0.543 

0 

0 

0.8 

27T 

256 

256 

256 

23.172 

23 

23 

47T 

256 

256 

256 

5.793 

5 

5 

87T 

256 

256 

256 

1.448 

1 

1 

107T 

225.456 

181 

181 

0.927 

0 

0 

167T 

55.043 

35 

35 

0.362 

0 

0 


The first differential approximation to the first order accurate driver operator (2.1) taken for the char- 
acteristic components is given by 

fda(l 1 ) = % — T (f + t)d m , 

2(1 + i 2 ) 2 


(2.18) 


d 1 


2\/T+1? 


W2 (u>2hj (t 2 + t j 


and the 6-accuracy penetration distance is 


(2.19) Si = -d! In(l-c). 

The first test we perform to validate the discrete half- space analysis from Sec. 2.1 and analytical 
expressions (2.17) and (2.19) for penetration distances of characteristic components. We calculate penetration 
distances of e = 0.01 accuracy for different values of t and u>2 for operators (1.4) and (2.1). Our aim is to 
compare the following distances from the boundary (measured in meshsizes) : (1) calculated by formulas 
(2.17) and (2.19) (<s£^ = 6j/y/l +t 2 /h, j = 1,2), (2) derived from the half-space analyses from Sec. 2.1 

) an d (3) computed in direct numerical simulations (<5^). In analytical calculations, the exact solution 
of the differential problem (1.1), (1.3) is always assumed to be e l ^~ tx ^ y \ In direct numerical simulations 
the exact solution has been chosen to be sin(cj 2 ( — tx + j/)). The numerical distance is considered to be m 
if the Lqo norm of the discretization error on (m + l)-th vertical line is greater than 0.01. The simulation 
grid is 257 x 257. If in analytical calculations the result exceeded 256 (the penetration distance covers all 
the domain) it has been set to 256. Table 2.1 contains the test results. The two obvious conclusions are 
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1. The discrete half-space analysis is actually precise. In all the tests, the results predicted by this 
analysis and obtained in real numerical calculations coincide. 

2. The inequality (2.16) provides a good estimate of the penetration distance, especially, for the first 
order operator or for small angles of attack (t < 0.6). Some deterioration in predicting the second 
order operator penetration distances for nearly diagonal alignment is explained by the fact that the 
penetration distance in case of 45° angle of attack (t = 1) is determined by the third order term 
which is not taken into account in calculating tfy Nevertheless, the estimate obtained from (2.16) 
seems to be reliable to predicting the key property: whether the penetration distance is comparable 
with (or larger than) the characteristic size of the domain. Moreover, it gives us a possibility to 
estimate the ratio between the penetration distances of the target and driver operators even on very 
coarse grids where one of the distances (or both of them) is shorter than one mesh size. 

2.3.2. Convergence to Within the Given Accuracy. Let function U(x,y) — e *(^i x+^y) b e the 
exact solution of (1.1) and (1.3). We can reformulate the discrete problem in the following way: we are 
looking for a discrete function Ui lt i 2 (an approximate solution to (1.4)) that possesses the total error (in the 
Loo norm) which is not greater than given c, i.c., 


max 

*1 j *2 


U{i\h,i2h) 


< €. 


In order to find such a solution, the discretization error of the target operator, of course, should satisfy 


( 2 . 20 ) 


max 

*l ,*2 


U(iih,i2h) 


< 


where is the exact solution of (1.4). This, in particular, implies that the target operator penetration 

distance 62 of the e - accuracy on the uniform grid with spacing h for incoming frequency U2 is larger than 
the characteristic size R of the domain (R = vl + t 2 in our case). The condition (2.20) actually defines the 
coarsest possible grid on which the desired e - accuracy can be achieved. The goal of converging within the 
discretization error, which is typical in FMG type algorithms, can be considered as a particular case where 
e is the target operator discretization error. 

The analysis from Sec. 2.1 provides us with an upper bound Q for frequencies of incoming Fourier modes 
LJ 2 satisfying (2.20). The same analysis predicts the number of defect-correction sweeps N sweeps required to 
achieve an approximation possessing the e-accuracy. Both a) and N sweeps depend, of course, on the mesh 
size h of the given grid, on the desired accuracy e, and on the non-alignment parameter t. 

Fig. 1 demonstrates the typical behavior of u> and N sweeps as functions of the mesh size h. In all the 
experiments we performed, the 1%-accuracy (c = 0.01) was picked on. The choice of angles of attack was 
representative. For simplicity, we always started from the initial approximation in the interior of the domain 
(fy >0) obtained from the solution of the driver operator. 

The results corroborate the conclusions made in Sec. 1.2 about the growth of the number of cycles 
required to achieve an approximation possessing the e- accuracy. For obviousness we added to the lower plot 
on Fig. 2.1 the solid line corresponding to the function fi -1 / 3 . 

The practical conclusions are the following. 

1. The number of DC1 iterations required to achieve a given accuracy grows as some (negative) power 
of h (approximately h~*). 

2. For any given (continuous) boundary conditions, the analysis is able to provide predictions of the 
grid required to solve this problem to a desired accuracy and how many defect-correction sweeps 
should be performed to achieve this accuracy. 
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Fig. 2.1. a j and N 3we e P s as functions of the meshsize h . 

3. If the flow direction is variable, one should select uj and Sweeps such that w satisfies (2.20) for all 
the possible angles of attack and Sweeps iterations provide the desired accuracy for characteristic 
components with the incoming frequency u for any possible t. 

4. A smart choice of the initial approximation (e.g., the initial approximation being interpolated from 
the coarse grid in the framework of an FMG solver) can sometimes reduce N sweeps but the qualita- 
tive behavior remains the same: the number of required defect-correction iterations grows as some 
negative power of h in passing to finer grids. 

There are several ways to change the algorithm in order to get N sweeps independent on the meshsize h : 

1. The first possibility which is studied in Sec. 3 is to apply a driver of the same approximation order 
as the target operator. This cure is actual only for second order operators since for higher-order 
scheme there are no stable upwind discretizations. 

2. The second way is to use a predictor- corrector technique for solving the target operator. This method 
suggests some marching along the flow direction. If this is possible (there are no recirculation zones), 
then this approach is, probably, one of the best allowing to achieve an optimal efficiency. 

3. The third method is a smart multigrid algorithm which employs the coarse-grid operators approxi- 
mating the characteristic component FDA of the target operator. 

The two latest approaches are the subjects of future papers. 

Before going to the defect-correction scheme with the second order driver we consider the asymptotic 
convergence of the first order driver scheme. 

2.3.3. Asymptotic Convergence. First of all we believe that the most important characteristic of the 
solver efficiency is the ability to fast solve the problem to the desired accuracy. The fact that the accuracy 
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t=0.6 




Fig. 2.2. DC1: residual convergence rate history. 


has already been reached can be established from the comparison of the solutions on grids with different 
meshsizes. In this view, delivering the residuals to the computer zero is less important. Other criteria (rather 
than the vanished residuals) should be used to decide on stopping iterations on the given grid and passing 
to a finer grid. However, the principal possibility to drive residuals to the computer zero is considered to be 
useful. Tests below demonstrate that DC1 iterations possess this property. 

It was observed by many researchers that the asymptotic convergence rate of the defect- correction solver 
for (1.4) is about 0.5 per iteration. (In [5], authors emphasized that the asymptotic convergence rate can 
somehow deteriorate for the central and the pure upwind target discretizations.) This asymptotic convergence 
rate can actually be further improved by a proper residual weighting ([1] suggests the weight 2/3 for the pure 
upwind operator). However, we already mentioned that this good convergence rate is only achieved after 
many sweeps with a much poorer convergence. Fig. 2.2 demonstrate the residual convergence history for 
different representative vales of the incoming frequency o »2 and the non-alignment parameter t. In the legend, 
the corresponding amplification factors (the asymptotic convergence rate (p(G r )) and the convergence rate 
bound (||G r || 2 )) calculated by the methodology explained in Sec. 2.2 are shown. In all the numerical tests 
performed for the homogeneous problem (1.4) {fi u i 2 = 0) on the uniform grid with meshsize h = 2 -6 , the 
iterations was stopped when the residual L 2 norm became less than 10 -10 . The results of the tests can be 
summarized as following. 

1. The asymptotic convergence rate is always good enough and the p(G r ) estimate is its accurate 
prediction. However, this good convergence is manifested only either on very fine grids or after a 
lot of iterations. 

2. The bound ||G r || 2 is not very sharp in the presented tests. The full space Fourier analysis (eliminating 
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the characteristic components from the consideration) gives a similar estimate (see [5]). It, probably, 
means that in order to observe this “worst” behavior one should test non-homogeneous problems. 

3. The number of sweeps required to get into the asymptotic convergence regime is roughly propor- 

tional to the ratio mm(R,d 2 )/d u where d x (2.18) and d 2 (2.15) are the first and the second order 
penetration distance parameters respectively. This number grows roughly as in passing to 

finer grids. 

4. For any given incoming frequency u> 2 , there is a grid with small enough meshsize h (d 1 = O(R)) on 
which the defect-correction iterations demonstrate the asymptotic convergence rate from the very 
beginning. 

3. Defect-Correction Method with Second Order Driver (DC2). The second order accurate 
upwind discretization is defined as 


L n u tui2 — 2^ 1 4- t 2 h (( 3uil -4u u-i,i2 + w h-2,* 2 ) 

*i = 1, 2, ... TV, i 2 = 1, 2, . . . TV; 
u o,i 2 = 9(i 2h), u -i,%2 = g'{hh). 

It is stable in marching and can serve as a driver for the defect-correction iterations. 

We tested the DC2 iterations for the same test-cases as in Sec. 2. In all the experiments on all the 
grids, the number of DC2 iterations required to get 1% accuracy did not exceed 3 (including the first sweep 
marching the driver operator (3.1) to obtain the initial approximation). The residual convergence rate history 
shown on Fig. 3.1 confirms the predicted efficiency of the defect-correction method with the second order 
driver and demonstrates the reliability of the residual convergence analysis introduced in Sec. 2.2. 

4, Adaptive Multigrid Algorithm (AMA). Any adaptive multigrid solver is usually required to 
make two important decisions. The first is to decide on stopping iterating on the given grid and proceeding 
to the next fine grid. The second issue is to realize that the required approximation is achieved and finish 
its work. 

Some prior information being incorporated into the solver can improve its performance a lot. For 
example, if the amplitudes of all the significant (for the given e-accuracy) Fourier components of the inflow 
boundary conditions can be estimated, then many calculations can be performed in advance. In particular, 
using the half- space analysis for the highest frequency essential Fourier component, one can find the optimal 
grid spacing h and the number of the defect- correct ion iterations on that grid required to achieve the desired 
accuracy. Then, a single-grid algorithm performing the necessary AT sweeps iterations seems to be the optimal 
solver. 

For general boundary conditions, this approach does not work. In case of general boundary conditions, 
however, all the decision should be done “automatically” employing only the data computed during the 
solution process. The algorithm we propose in this section is based on the full multigrid methodology where 
the comparison between solutions on different grids becomes a criterion for stopping further refinement. 
The adaptive multigrid algorithm solving the problem to the e-accuracy can be recursively defined in the 
following 5 steps; 
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1 = 0.2 



Fig. 3.1. DC2: residual convergence rate history. 

1. Let u 2h be the solution on the grid with meshsize 2 h. On the grid with meshsize h for the given e- 
accuracy, the highest possible incoming frequency u > for which this accuracy is achievable is calculated 
from the half-space analysis. The required number of the sweeps AT swceps is defined as well. 

2. The initial approximation on the interior of the h-grid is obtained by (bilinear) solution interpolation 
from the 2h-grid. 

3. One defect-correction iteration is performed and an approximate solution u h is obtained. 

4. If the L 0 o norm of the difference between u h and u 2h is less than e then u h is the final solution. 
Solver finishes its work. 

Otherwise, additional defect-correction sweeps are performed. The iterations on the h- grid are 
stopped when either the L ^ norm of the difference between the two successive approximations 
becomes less than t or the total number of sweeps (including the first one performed in the Step 3) 
reaches the corresponding h-grid N sweeps . In fact, as one can see from the numerical tests below, 
the latter tolerance was reached only on relatively coarse grids, where the target discretization has 
no accuracy at all. The last approximation is considered as the /i-grid solution u h . 

5. h replaced with h/2. Go to Step 1 for the next fine grid. 

Table 4 collects u> and iV 8weeps values on grids with different meshsizes h for e — 0.01. In this table, 
u) is the highest frequency resolved in the target second order operator on the h- grid to the e- accuracy for 
any angle of attack (0 < t < 1) and iV gweep8 is the number of DC1 sweeps ensures the convergence within 
this accuracy. We found from the analysis and checked in numerical tests that on grids with h > 2~ 7 the 
maximal cross- characteristic interaction defining u> is observed for the diagonal alignment case (t = 1) while 
on the finer grids the strongest interaction occurs at t « 0.6. The maximal AT sweeps is always found at the 
diagonal alignment. 

If the driver is second order accurate then on all the grids the number of required DC2 sweeps is bounded 

to N sweeps — 3. 

The Figs. 4.1 and 4.2 demonstrate the performance of AMA based on DC1 and DC2 defect-correction 
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Table 4.1 

Limits for iterating DC 1 scheme on different grids. 


h 

2~ 4 

2~ 5 

2-6 

2~ 7 

2-8 

2 -9 

2-1° 

2 -n 

2-12 

2~ 13 

2~ 14 

u 

1.17T 

1.97T 

3.2 7T 

5.47T 

8.6tt 

13.7tt 

21.87T 

34.7tt 

55.1t r 

87.6tt 

139.17T 

N 

1 ’ sweeps 

5 

6 

7 

9 

9 

10 

12 

14 

16 

18 

21 


t = 0.8 



Fig. 4.1. DCl: Adaptive multigrid algorithm (DCl-AMA). 


iterations. The input data were chosen so that the function U(: r, y) = sin^u^y- tx)J is the exact solution of 
(1.1), (1.3). The non-alignment parameter was set to t = 0.8. We tested different frequencies u> 2 providing a 
smooth solution to the problem with y-periodic boundary conditions for which the 1%-accurate solution can 
be obtained on a grid with h > 2~ 9 . The vertical coordinate on the figures marks (logarithm of) the total 
error and vertical lines separate the results corresponding to calculations on different grids. The first value 
on each grid (except the coarsest grid) is the value of the total error after the solution interpolation from the 
previous coarse grid. All the next values indicate the total error after the corresponding defect-correction 
iterations. The adaptive algorithm proved to be quite efficient requiring one extra level iteration at most to 
ascertain that the 1%-accuracy is already achieved. This is a very reasonable cost equal to the cost of few 
additional sweeps on the coarsest possible grid where the 1%-accuracy could be reached. 

The small number of iterations performed by DCl-AMA on the fine grids does not disprove the claim 
that the number of required iterations might grow on the fine grids. The accuracy considerations make sense 
only for differentiable solution. The three solution conditions (the vertical periodicity, the 1%-accuracy on 
a grid with h > 2 9 and the differentiability) together leave us just a few allowed components (only those 
with uj 2 — 27 r/c, k = 0, 1, 2, 3, 4, 5, 6). Therefore, we have no chance to approach on the tested grids the 
component which realizes this predicted behavior. Of course, on finer grids the expected behavior is much 
likely to be manifested. 

The efficiency of DC2-AMA seems to be nearly optimal. 
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t = 0.8 



Frc. 4.2. DC2: Adaptive multigrid algorithm (DC2-AMA). 

5. Conclusions. To summarize the practical results of the research reported in this paper we make 
the following conclusions: 

1. The efficiency of the defect-correction method with the first order driver is grid dependent. The 
number of iterations required to reach the desired accuracy (or the asymptotic convergence rate 
regime) might grow on the fine grids (roughly as h~ 5). 

2. Using the second order driver in the defect-correction iterations eliminate this dependence resulting 
in a very efficient solver. We are aware that in many cases the choice of the first order driver 
is dictated by external reasons. For example, in solving discretized multidimensional hyperbolic 
systems of equations where downstream marching is impossible, first order schemes are considered 
to be much easier to solve than higher-order schemes. Nevertheless, we believe that even in such 
cases the opportunity of employing the second order driver should be carefully studied. 

3. Any robust solver using defect-correction iterations should adopt the adaptive multigrid approach. 
Several further simplifications can be suggested on this way. For example, if the boundary conditions 
and/or the problem geometry prove to be too complicated, so that it seems hard to estimate N 8yfeeps 
on each grid, then the second criterion in Step 4 for stopping iterations on a given grid (performing 
all AT S W eeps iterations) can be omitted. In many cases this results in some additional work on coarse 
grids but does not affect the total work count. 

4. The discrete half-space analysis described in Sec. 2.1 is an accurate and very efficient tool for 
predicting actual solution behavior. 

5. The defect- correction iterations converge residuals to the computer zero. The asymptotic conver- 
gence rate and an upper (but not the sharpest) bound can be calculated by means of the matrix 
analysis reported in Sec. 2.2 
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order, then the initial convergence rates may be very slow. The number of iterations required to get into the 
asymptotic convergence regime might grow on fine grids as a negative power of h. In the case of a second order 
target operator and a first order driver operator, this number of iterations is roughly proportional to A - */®. (3) If 
both the operators have the second approximation order, the defect-correction solver demonstrates the asymptotic 
convergence rate after three iterations at most. The same three iterations are required to converge algebraic error 
below the truncation error level. A novel comprehensive half-space Fourier mode analysis (which, by the way, can 
take into account the influence of discretized outflow boundary conditions as well) for the defect-correction method 
is developed. This analysis explains many phenomena observed in solving non-elliptic equations and provides a close 
prediction of the actual solution behavior. It predicts the convergence rate for each iteration and the asymptotic 
convergence rate. As a result of this analysis, a new very efficient adaptive multigrid algorithm solving the discrete 
problem to within a given accuracy is proposed. Numerical simulations confirm the accuracy of the analysis and the 
efficiency of the proposed algorithm. The results of the numerical tests are reported. 
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